Deducing Local Rules for Solving Global 
Tasks with Random Boolean Networks 



Bertrand Mesot ^ Christof Teuscher '^'^ 

^IDIAP Research Institute, Rue de Simplon 4, CH-1920, Martigny, Switzerland 

^Los Alamos National Laboratory, Advanced Computing Laboratory, CCS-1, 
MS-B287, Los Alamos, NM 87545, USA 



Abstract 

It has been shown that uniform as well as non-uniform cellular automata (CA) can 
be evolved to perform certain computational tasks. Random Boolean networks are 
a generalization of two-state cellular automata, where the interconnection topology 
and the cell's rules are specified at random. 

Here we present a novel analytical approach to find the local rules of random 
Boolean networks (RBNs) to solve the global density classification and the synchro- 
nization task from any initial configuration. We quantitatively and qualitatively 
compare our results with previously published work on cellular automata and show 
that randomly interconnected automata are computationally more efficient in solv- 
ing these two global tasks. Our approach also provides convergence and quality es- 
timates and allows the networks to be randomly rewired during operation, without 
affecting the global performance. Finally, we show that RBNs outperform small- 
world topologies on the density classification task and that they perform equally 
well on the synchronization task. 

Our novel approach and the results may have applications in designing robust 
complex networks and locally interacting distributed computing systems for solving 
global tasks. 
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1 Introduction 



Cellular automata (CA) [51,55] were originally conceived by Ulam and von 
Neumann [53] in the 1940s to provide a formal framework for investigating 
the behavior of complex, extended systems. CAs are dynamical systems in 
which space and time are discrete. A cellular automaton usually consists of a 
D-dimensional regular lattice of N lattice sites, commonly called cells. Each 
cell can be in one of a finite number of S possible states and further consists of 
a transition function F (also called rule), which maps the neighboring states 
to the set of cell states. CAs are called uniform if all cells contain the same 
rule, otherwise they are non-uniform. Each cell takes as input the states of the 
cells within some finite local neighborhood. By convention, a cell is considered 
to be a member of its own neighborhood. In the standard one-dimensional 
CA model, for example, each cell is connected to the r (r stands for radius) 
immediate local neighbors on either side as to itself, thus each cell is connected 
to 2r + 1 neighbors. 

Random Boolean networks (RBNs), on the other hand, form a more general 
class of discrete dynamical systems, in which two-state cellular automata are 
a special subclass. Random networks and RBNs were investigated in the past 
forty years by many a researcher (see for example [2,4,7,36]), but became 
popular mainly due to the contributions of Stuart Kauffman [24-26] and oth- 
ers. 

In its simplest form, a RBN is composed of nodes (sometimes also called 
elements or cells) where each node can be in one of two possible states (0 or 
1) and receives inputs from K randomly chosen other nodes (self-connections 
are allowed). Note that K can refer to the exact or to the average number 
of connections between the nodes. In the geneticist's term, for example, K 
measures the richness of epistatic interactions among the components of a 
system [26]. 

The deterministic behavior of each node is specified by one out of the 2^^ 
Boolean functions, specifying its next value for each of the 2^^ combinations 
with K inputs. Like CAs, RBNs can be non-uniform (i.e., each node can 
potentially have a different rule) or uniform, although they are non-uniform 
in the majority of the cases. The network's nodes are usually updated syn- 
chronously, although multiple asynchronous updating schemes exist (see [20] 
for an overview). 

Synchronous random Boolean networks as introduced by Kauffman are com- 
monly called NK networks or models. Figure 1 shows a possible NK random 
Boolean network representation (A^ = 8, K = 3). The model is very similar 
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to the well-studied class of models which arises in statistical physics, called 
spin-glasses [18]. Spin-glasses are disordered magnetic materials in which the 
orientation of nearby magnetic dipoles may be either parallel or antiparallel 
in space. 




rule 43 



Fig. 1. Illustration of a random Boolean network with N = 8 nodes and K = 3 inputs 
per node (self-connections are allowed). The node rules are commonly represented 
by lookup-tables (LUTs), which associate a 1-bit output (the node's future state) to 
each possible K-hit input configuration. The table's out-column is commonly called 
the rule of the node. 

More formally speaking, let s* be the state-vector of the network's nodes at 
time t. The neighborhood Il{xi) of a node Xi is defined as 

U{xi) C Pk{{xi, . . .,xn}), 

where Pk{X) denotes the set of all subsets of size K from X (self-connections 
are allowed). The neighborhood size is given hy K = |n(xj)|. The most com- 
mon ways to choose the neighborhood are the following: 

• random neighborhood: choose the K neighbors randomly from the set 
{xi, . . . , xn} \ {xi\ (without self-connections) or from the set {xi, . . . , xjq} 
(self- connect ions allowed) 

• adjacent neighborhood: K neighbors are randomly chosen in the "immedi- 
ate" neighborhood. 

A(7V, K) is sometimes used to denote a iVK model with adjacent neighborhood, 
N{N, K) for the random neighborhood. We shall not go into more details of the 
adjacent neighborhood here, as we will only consider purely random networks 
with a uniform probability to connect two nodes together in this article. 

Figure 1 shows the interaction graph G{y^ E) where V = {xi, . . . , x^} corre- 
sponds to the set of nodes, and {xi,Xj} G if and only if Xi is connected to 
Xj. The node rules are commonly represented by lookup-tables (LUTs), which 
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associate a 1-bit output (the node's future state) to each possible K-hit input 
configuration. The table's out-column is commonly called the rule of the node. 



Whereas the homogeneous connectivity and the usually uniform rules of a 
CA restrict the automaton's parameter space, RBNs have a vastly greater 
space, which adds an additional challenge for both evolutionary and analytical 
methods, either to analyze the networks' behavior or to find networks able to 
perform a given task. 



The main contribution of this paper consists in a novel analytical method 
which allows to determine the local rules of a random Boolean network for 
two global and well-known problems for cellular automata, (1) the synchro- 
nization and (2) the density classification task. Whereas previous work was 
purely quantitative and focused on the co-evolution of cellular automata to 
find suitable rules and architectures [45], our qualitative approach goes a step 
further and also provides convergence and quality estimates. We systematically 
compare the results with previous work on standard and non-standard cellu- 
lar automata architectures, including small-world networks, and show that 
randomly interconnected automata perform better than a standard cellular 
automata with a regular interconnection architecture. Our method also allows 
the networks to be randomly rewired during operation without affecting the 
global performance. The approach and results may have applications in de- 
signing robust complex networks and locally interacting distributed computing 
systems for solving global tasks. 



In Section 2 we describe the two global tasks used in this paper, the density 
classification task and the synchronization task. The main challenge in solving 
such global tasks with CAs and RBNs consists in finding the node's local rules 
that will result in the desired global automata behavior. Note that no gen- 
erally applicable method for this exists. In Section 3 we derive the recursive 
equation which gives the probability of a node to be in state 1 at the next 
time step. This equation is then used in Section 4 to determine the node's rules 
for both tasks. Section 5 examines the network's performance as a function 
of N and K whereas Section 5.3 compares the results with previously pub- 
lished performance results of cellular automata for the same tasks. In order to 
measure the rule's capacity to solve the tasks, we introduce an entropy-based 
performance measure in Section 5.4. The same measure also allows us to pre- 
dict how quickly the two tasks can be solved. In Section 6 we briefiy discuss 
non-uniform random Boolean networks and perform experiments with small- 
world network topologies. Our findings, their possible future applications, and 
future work are discussed in Section 7. 
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2 Task Definitions 



Two commonly used applications for CAs are the density and the synchroniza- 
tion task. Both of these "global" tasks are mostly trivial to solve if one has a 
global view on the system (i.e., if one has access to the state of all nodes at the 
same time), but are non-trivial to solve for CA and RBNs, mainly because of 
the locality and the limited number of their interconnections, the simplicity of 
the basic cells (i.e., the basic components), and because there is no generally 
valid way to determine the cells's rules for such a massive parallel system. 
Unlike in the standard approach to parallel computation, in which a given 
problem is split into independent sub-problems, CAs and RBNs have to solve 
problems by the bottom-up approach, where the global and usually complex 
behavior arises from nonlinear, spatially extended, and local interactions [30]. 
This difficulty has naturally also limited CAs applications. 

In the following, the density classification and the synchronization task shall 
be briefiy described. 



2.1 The Density Classification Task 



The density classification task for CAs must decide whether or not the initial 
configuration of the automaton contains more than 50% of Is. In this context, 
the term "configuration" refers to an assignment of the states or 1 to each 
cell of the CA (i.e., there are 2^ possible initial configurations). The desired 
behavior of the automaton is to have all of its cells set to 1 if the initial 
density of Is exceeded 1/2, and all otherwise. The density task for CAs can 
straightforwardly be applied to RBNs. The special case of having an equal 
number of Is and Os in the network is commonly avoided by using an odd 
number of nodes. 

The density task was studied among others by Mitchell et al. [29-31], Das et 
al. [15,16], Sipper [38-41], Sipper et al. [43,45,46], Capcarrere [14], Capcarrere 
et al. [10-12] and Tomassini et al. [52] in various forms, including non-uniform 
CAs, asynchronous CAs, and non-standard architectures. It should also be 
mentioned that no uniform two-state CA exists, which perfectly solves the 
density classification task for all initial configurations [27]. 

Let n\ be the number of nodes being in state 1 at time t. The density classi- 
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Fig. 2. Three possible maps which satisfy Equation 1. The dash-dotted curve rep- 
resents the one that wiU take the highest number of iterations before settUng down 
in a fixed point. 

fication task can then be formally described by the following equation: 



3tlim S.t. t ^ tli) 



if n° < N/2 

iV if > N/2 
N 



I 2 



otherwise 



Let us define p\ = n\/N as the fraction of nodes that are in state 1 at time 
t. We will sometimes use p\ instead of n\ when the networks are sufficiently 
large (i.e., N — >■ +oo), since p\ may then be interpreted as the probability 
of the node to be in state 1 at time t. Moreover, since we only consider a 
synchronous updating scheme, p*^^ will solely depend on p\. The network's 
behavior can thus simply be described by a map of the form p^^^ = Q{p\), 
where Q : M i — > M. Each of the three solutions of the density classification 
task will then be represented by a fixed point in this map. Figure 2 shows 
three possible functions Q that satisfy Equation 1. Note that the fixed point 
inpl = 1/2 is unstable, which means that a small perturbation of the network's 
state, such as for example one state change among the nodes, will drive the 
system towards the wrong fixed point, i.e., towards the wrong solution. This 
unstable fixed point can fortunately be avoided by using an odd number of 
nodes. 
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2.2 The Synchronization Task 



The one- dimensional synchronization task (also called firefly task) for syn- 
chronous CA was introduced by Das et al. [15] and studied among others by 
Hordijk [21] and Sipper [42]. In this task, the two-state one- dimensional or 
two-dimensional CA, given any initial configuration, must reach a final con- 
figuration within M time steps, that oscillates between all Os and all Is on 
successive time steps. The whole automaton is then globally synchronized. 

As with the density classification task, synchronization is a non-trivial global 
computation task for a small-radius CA where all cells must coordinate their 
behavior with all the other cells while having only a very local view of their 
neighborhood. The two-state machine of each cell obviously also prevents any 
counting, which would make the task much less difficult. Furthermore, in the 
non-uniform case, where each cell can have a different rule, there is an immedi- 
ate solution consisting of a unique "master" rule that alternates between and 
1 without looking at his neighbors, and all other cells being its "slave", i.e., 
alternating according to their right (or upper, in the two-dimensional case) 
neighbor state only. However, Sipper [42] used non-uniform CAs to find per- 
fect synchronizing CAs by means of evolution only. It appeared that this very 
particular solution was never found by evolution and, in fact, the "master" 
or "blind" rule 10101010 (rule 170) in one dimension, was never part of the 
evolved solutions. This is simply due to the fact that this rule has to be unique 
for the solution to be perfect, which is contradictory to the natural tendency 
of the evolutionary algorithm used, as was demonstrated by Capcarrere [12]. 

The two-dimensional version of the task is identical to the one- dimensional 
version, except that the necessary number of time steps granted to synchronize 
is not anymore in the order of A^, the number of cells in the automaton, but 
in the order of n -|- m, where n and m are the size of each side of the standard 
CA. Therefore, the speed of synchronization is much faster compared to the 
one-dimensional version. 

Extending the two-dimensional synchronization task for CAs to RBNs is again 
straightforward and solutions found be means of co-evolutionary algorithms 
were presented by Teuscher and Capcarrere in [50]. The firefiy synchronization 
task was also successfully implemented in actual evolutionary hardware, for 
both the one [44] and the two-dimensional [50] version. 

Similarly to the density classification task, the synchronization task may also 
be formally described by a condition on the number of Is in the network's 
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Fig. 3. Three possible maps which satisfy Equation 2. The dash-dotted curve requires 
the highest number of iterations before oscihating between all pi = and pi = 1. 



Figure 3 shows three possible maps which satisfy Equation 2. Contrary to the 
density classification task, fixed points must be avoided since oscillations are 
required. 



3 Prom Global Configurations to Local Interactions 

The challenge in solving the density classification and synchronization task 
with CAs and RBNs consists in finding the node's local rules that will result 
in a global automaton behavior which satisfies Equations 1 or 2 for all possible 
initial configurations. The difficulty comes from the fact that each node has 
only access to K randomly chosen bits of the entire RBN's state vector s* at 
time t. It is therefore crucial to know how the global network state is seen by 
the local nodes in order to understand the network's dynamics. 

Let Wl be a random variable representing the number of Is in a node's neigh- 
borhood at time t and let and = N{/N be two random variables which 
represent the number of Is and the percentage of Is in the network's state s* at 



state: 




(2) 
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time t respectively. obviously depends on P/, thus, if we suppose that the 
nodes' neighbors are randomly chosen from a uniform distribution, then the 
probability of having d neighbors in state 1 at time t, knowing the percentage 
of Is (i.e., p\) in the global network state s* of a RBN, is therefore given by: 

P(wi = d I p* = p\) = (^^^ {p\r (1 - p\)^-'. (3) 

The distribution of Wl knowing Pi thus follows a binomial law B{K,p\) and 
the expected number of Is in a node's neighborhood is therefore Kp\, which 
implies that the average number of Is in the input of a node is the same as the 
average number of Is in the global network state. This property is actually 
the key that will allow us to reduce the complex behavior of TV nodes to the 
stochastic behavior of a single node which behaves on the average like the 
entire network. 

Let s* be the state of a RBN at time t, composed of nodes, each being 
randomly connected to K other nodes (self-connections are allowed). Further- 
more, let us assume that N is sufficiently large to consider p\ = n\/N as the 
probability of a node to be in state 1 at time t. We will later see in our ex- 
periments that this approximation only affects performance. Knowing s*, the 
probability that the future state of node z is 1 is then given by: 

P(5*+i = 1 1 S* = s*) = P(5*+i = 1 1 = j) P(<f * = J I S* = s*), (4) 

j=0 

where S'*"'"^ and $* are two random variables that represent the state of node i 
at time t + 1 and its input at time t respectively. This is a probabilistic descrip- 
tion of the node's Boolean transfer function (i.e., its rule), where P($* | S*) 
and P{Sl~^^ I $•) are the input and output distributions respectively. If fi{j) is 
the output of node i for an input value j, then the previous equation becomes: 

P(5*+i = 1 1 S* = s*) = m P{¥, = J I S* = s*), (5) 

j=0 

since the probability of a node to be in state 1 at time t + 1 is 1 if fi{j) = 1, 
and otherwise. By using the variables and P* as introduced above, this 
latter equation can further be reduced to: 

PiSl-"' = 1 1 Pi = P\) = E c, P{Wl = J I Pi = p\), (6) 

i=o 

where j is the number of Is in the input of node i and Cj is the probability that 
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the future state (the output) of node 2 is 1, provided that the input contains 
j times a 1. 

Putting Equation 3 into Equation 6 gives us: 

p{sr = 1 1 Pi = p\) = E c, f f ) {p\y (i - pi)''"^' , (7) 

,=0 \J/ 



t+1 



which holds for any node i. Since P{Sl~^^ = 1 | P* = p*) is equivalent to 
the probability for a node to be in state 1 at time t + 1, we can now express 
Pi"^^ as a function of p\ by means of the following recursive equation: 



K 



ri'^ = E«i(pi) (i 

j=0 



Pi 



K-j 



with a, 




where the parameters aj represent the number of the rule's inputs that contain 
j times a 1 and that produce a 1 as output. Note that, if we represent a 
node's rule by a lookup-table, then aj corresponds to the number of entries 
that contain j times a 1 and that have a 1 as output. For example, the two-bit 
XOR rule (i.e., the output is 1 if and only if the two inputs are different) 
corresponds to oq = 0, ai = 2 and a2 = since the output is 1 if there is only 
one 1 in the input. On the other hand, the two-bit NAND rule corresponds 
to: ao = 1, fli = 2 and 02 = 0. 

Recently, Matache and Heidel [28] used a formula to describe the probability 
of finding a node in state 1 at time t to analyze deterministic chaos in random 
Boolean networks. The basic idea is similar to our approach, but their work, 
which extends the model studied by Andrecut and Ali [6], is limited to rule 
126 only, whereas our approach is valid for any rules. 



4 Finding the Node's Rules 



4-1 Rules for the Density Classification Task 



As seen in Section 2.1, solving the density task can be reduced to finding maps 
Pi^^ = Qk{p\) which satisfy the constraints that were informally represented 
in Figure 2. More formally speaking, these constraints may be defined by the 
following set of equations: 

(CI) Qx(o) = o, gx(i) = i, gi^(i/2) = 1/2, g^(i/2) = o 
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(C2) Q';^{x) > Vx G [0, l/2[ and Q'^{x) < Vx G ]l/2, 1]. 



However, according to Equation 8, we have: 



K 



Qk{x) = ^ ttj x^ (1 — x)^ ^ with ttj G 



0, 



.J J 



(9) 



Applying constraints CI and C2 to Equation 9 and solving it as a function 
of Qj provides us now with an elegant means to directly determine the node's 
rules for the density task. In the following sections, different possible solu- 
tions for small connectivity parameters K shall be presented. We will see that 
the simple majority rule presents a common solution for all networks with a 
connectivity K > 2. 



4.1.1 K = 2: No Solution 
For K = 2, Equation 9 becomes: 

Q2{x) = ao(l — xY + aix{l — x) + a2X^ . 

After applying constraints CI we obtain the following set of values: 



Q2(0) = ^ ao = 
Q^il) = 1 ^ a2 = 1 

Q2(l/2) = ^(ai + l) = ^ ^ ai = l 

This implies Q2{x) = x and therefore Q2{x) = for all x, which does not 
satisfy constraints C2. Hence, for = 2 no set of rules is able to perfectly 
solve the density classification task for all initial configurations. 



4.1.2 K = 3: A First Solution 
For K = 3, Equation 9 becomes: 

Qsix) = ao(l - xj^ + aix{l - x)"^ + a2X^{l - x) + a^x^. 



In order to satisfy constraints CI, we must impose ao = 0, 03 = 1 and 02 = 
3 — Oi. Qsi^x) and Q'^^x) then become: 
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Q3{x) = 2{ai-I)x^ + 3{1 
Q'^{x) = 12{ai - l)x + 6(l 



ai)x^ + aix 
ai). 



After applying constraints C2, 



we obtain: 



Q'^{x) > 12 (ai 



l)x > 6(ai - 1) =^ a; < - iif ai = 0. 



One can easily verify that oi = is the sole solution. The corresponding 
polynomial is therefore: 



This solution simply corresponds to the well-known majority rule: its output 
is 1 if and only if there are more Is than Os in the inputs, and otherwise. As 
one can verify, since oq = ai =0, the only input with all bits set to and the 



j = 3 inputs with one bit set to 1 result in a as output. Conversely, for 
all inputs containing more than one bit set to 1, the output becomes 1, since 
02 = 3 and 03 = 1. 

4.1.3 K = 4: Many Rules 

For K = 4 only one polynomial satisfies both sets of constraints. It corresponds 
to: ao = 0, ai = 0, 02 = 3, 03 = 4 and 04 = 1. Since 02 = 3 and according to 
Equation 9, of the (^^^ = 6 inputs with two bits set to 1, only half of them will 
set the output to 1. Since there are 20 possibilities of setting 3 out of 6 outputs 
to 1, there are 20 (out of 2^^ = 65536) equivalent rules which will successfully 
solve the density classification task. These rules are represented in Table 1. As 
one can see, all rules are in a certain sense majority rules, but applied to an 
even number of connections, reason why many possibilities exist. 

4-1-4 K = 5: Many Polynomials 

By the same procedure we find the values of oj when K = 5: 

ao = 0, ai = 0, 02 G [0, 5], 03 = 10 — 02, 04 = 5 and a^ = 1. 

Since 02 is the only free parameter, 6 polynomials satisfy conditions CI and 
C2, each corresponding to a value of 02- The majority rule is the only possible 
rule when 02 = 0, whereas in the other cases, each polynomial represents many 
rules. For example, when 02 = 3, (3°) (7°) = 14,400 equivalent rules exist. 



Q3{x) = 3x^(1 - x) +x^. 



(10) 
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1110111010100000 
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Table 1 

The 20 (out of 2^^ = 65536) rules for K = A random Boolean networks that suc- 
cessfully solve the density classification task for all initial configurations. 

As one can see, representing the rules by means of a polynomial bears the 
advantage that rules with similar dynamics are grouped together by Equation 
9. They may thus more easily be handled and analyzed on a meta-level. 



4-2 Rules for the Synchronization Task 

The procedure to find the rules for the synchronization task is similar to the 
density classification task. Again, we are looking for maps p*^^ = Qk{Pi) 
similar to those shown on Figure 3. More formally, the constraints we have to 
impose on Qk{x) are firstly 

Qk{0) = 1 and Qk{1)=0, 

which will force the network to oscillate between all Os and all Is on con- 
secutive time steps. And secondly, in order to make sure that the network's 
state converges towards = or = 1, Q'k{x) needs to be monotone and 
Q'^(x) ^Ofor allxe [0,1]. 

Two functions satisfying these constraints are for example: 

K-l 

QK{x) = {l-x f and gi^(x) = J] (1 - 

i=0 
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These functions are symmetrical and therefore perform equally well. We will 
call them 7;^-rules. Their behavior is very simple: the lookup-table associated 
with this rule outputs a 1 (respectively 0) if and only if all input bits are 
(respectively 1). Examples are rule 73 = 1 and rule 73 = 128 for K = ?> net- 
works. Both rules are basically fully equivalent and solve the synchronization 
task in the same way. 



5 Experiments and Results 

In this section we will examine the performance of RBNs as a function of 
the number of nodes and the number of connections per node K. Since 
the number of initial configurations exponentially grows with A^, testing them 
all becomes rapidly computationally intractable. Algorithm 1 illustrates the 
procedure we used for all experiments in this section. Since the nodes are ran- 
domly connected, we average the simulation results on a certain number of 
randomly generated networks. In addition, each initial configuration is gener- 
ated according to a biased distribution, i.e., the number of Is it contains is 
randomly chosen from a uniform distribution and the Is are then assigned to 
randomly chosen nodes. This will force the network to be tested uniformly on 
all possible values of p\ and not only around p\ = 1/2, which is the case if each 
configuration is chosen according to an unbiased distribution, i.e., each node 
has a uniform probability to be in state 1 or 0. Finally, in order to compare 

^min = 9 (minimal number of nodes) 
A^max = 201 (maximal number of nodes) 
rep = 5000 (number of initial configurations) 
nnet = 500 (number of networks) 

for n e [A^min, A^max] do 

for i G [1, nnet] do 

Generate a random network composed of n nodes, 
for r G [1, rep] do 

Generate an initial configuration. 

Test the configuration on the network, 
end for 

Average the results over rep. 
end for 

Average the results over nnet. 
end for 

Algorithm 1. Performance Evaluation Algorithm 

RBNs with CAs, we shall only test networks with an odd number of nodes. 
This avoids the unstable fixed point in = 1/2, which will be misclassified 
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Fig. 4. Density classification task performance as a function of and K. (1) solid 
line: K = 3, dashed line: = 4, (2) bottom-up, K = 5,7, 9, (3) bottom-up, 
K = 8,6. 

most of the time anyway, since, as explained in Section 2.1, a small perturba- 
tion, such as changing the state of one node, might be sufficient to drive the 
system towards the wrong solution. 



5.1 The Density Classification Task 



Figure 4 shows the network performance for K G {3,4,5,6,7,8,9}. All net- 
works were uniform and simply used the majority rule, which, as shown in 
Section 4.1, is the only common solution when K ^ 3. 

The results show that the performance increases with increasing A^. We will 
later show that it tends to 100% when — > +oo. Moreover, for odd values of 
K, RBNs perform better than for even values, and the larger K, the better 
the performance. This finding confirms that the majority rule works better 
without ambiguous situations, i.e., when a node does not receive an equal 
number of Is and Os. 

Figure 5 shows the average number of iterations a RBN requires to converge to 
a solution. Networks that need a small number of iterations to solve the density 
classification task are said to perform better. One can see that the number 
of iterations required increases with the network's size A^. This intuitively 
makes sense as a larger network needs more time to solve the task because 
the information must be propagated among the nodes. However, it is worth 
mentioning that the number of iterations required does not increase linearly 
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Fig. 5. Average number of iterations required to solve the density classification task 
as a function of N and K. (1) solid line: K = 3, dashed line: K = 4, (2) bottom-up, 
K = 9,7, (3) bottom-up, K = 5,8, 6. 

with A^, as one might expect, but rather follows a logarithmic-like law. This 
can be explained by the fact that random networks allow for long-distance 
connections, which help to propagate information faster. 

The performance increase with increasing can be explained by means of 
Equation 8: although the polynomials are continuous functions, we evaluate 
them at discrete time steps, i.e., p\ = n\/N can only take a finite number 
of discrete values. For example, let us suppose a network with = 10 and 
K = ?i. From a theoretical point of view, if p\ = 6/10, then p\^^ = 0.648, 
however, since p*"*"^ G {V-^ N ^ [0; -^]} the exact value taken by p^^^ will then 
depend on the actual network wiring. Performance will therefore get better 
with increasing and tends to perfection when A^ +oo. In order to further 
illustrate this, we generated Figures 6, 7, 8, and 9, that include both the 
theoretical and the experimental results. 

From these four figures we can see that the bigger A^, the better the theoretical 
curves are matched. As expected, we observe large variances for small values 
of A^, which implies that the network may converge toward the wrong solution 
near the fixed point p\ = 1/2, since p\^^ may cross the identity function 
Qk{,x) = X. This explains why the density classification task is solved less 
successfully for small A^. 

Moreover, the performance for odd values of K is better because the derivative 
at the fixed point p\ = 1/2 tends to be infinite when K ^ N . Errors in that 
point are then less likely to make p^"*"^ cross the identity function Qk{x) = x, 
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Fig. 6. Density classification task map p^^^ = Qk{Pi) for = 53 and K = 3. Mean 
(diamond) simulated and theoretical curve (dashed line). The identity function is 
also shown (dashed line). The bars indicate the maximum and minimum values 
observed. 




Fig. 7. Density classification task map p^"*"^ = Qk{Pi) for = 53 and K = 5. 
Mean (diamond) simulated values and theoretical curve (dashed line). The identity 
function is also shown (dashed line). The bars indicate the maximum and minimum 
values observed. 
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Fig. 8. Density classification task map p^^^ = Qk{p\) for N = 501 and K = ^. 
Mean (diamond) simulated values and theoretical curve (dashed line). The identity 
function is also shown (dashed line). The bars indicate the maximum and minimum 
values observed. 




0.2 0.4 0.6 0.8 1 



Fig. 9. Density classification task map p'l^^ = Qk{Pi) for N = 501 and K = b. 
Mean (diamond) simulated values and theoretical curve (dashed line). The identity 
function is also shown (dashed line). The bars indicate the maximum and minimum 
values observed. 
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Fig. 10. Synchronization task performance as a function of N and K. Bottom-up, 
K = 2,3,4,5,6,7,8,9. 

and therefore to drive the network towards the wrong solution. This can easily 
be seen if Figures 8 and 9: for K = 3, there are at least four points in the 
neighborhood of pi = 1/2 with an extrema that crosses the identity function, 
while there are only two for K = 5. Finally, note that the density classification 
task is trivially solved in one time step by the majority rule when K = N. 



5.2 The Synchronization Task 



In order to evaluate the performance of RBNs on the synchronization task, we 
used the following set of parameters (see Algorithm 1): A^'min = 10, A^ax = 200, 
rep = 10, 000 and nnet = 200. However, the initial configurations were still 
chosen according to a biased distribution as defined above. 

Figure 10 shows the performance of the synchronization task, whereas Figure 
11 shows the number of iterations necessary to synchronize the automaton. 
The network was uniform and all nodes used the 7/f-rule as described in 
Section 4.2. One can see that, similarly to the density classification task, per- 
formance increases with increasing A^ and K and that it tends to perfection 
(i.e., 100%) when A^ — > +oo or K ^ N. However, perfect results are already 
obtained for A^ ?a 150 and K > 3. Note that if A' = A^, each node possesses a 
full view on the entire network state and the task therefore becomes trivial. 



19 




20 40 60 80 100 120 140 160 180 200 
number of nodes 



Fig. 11. Average number of iterations required before synchronization occurs as a 
function of N and K. Top-down, K = 2, 3, 4, 5, 6, 7, 8, 9. 

5.3 RBNs versus CAs 

The performance of CAs for both the density classification and the synchro- 
nization task were investigated in various pubhcations (see Section 2.1 and 
2.2). We shall compare these results in the present section to those obtained 
with RBNs and will especially focus on how the information is processed in 
CAs and RBNs. 

In [29], Mitchell et al. presented four different CA rules for the density classi- 
fication task, namely ^gkl, (pmap (pcxp, and 0par, which all used a CA neigh- 
borhood of radius r = 3, i.e., each cell is connected to its 2r nearest neighbors 
and to itself. 0maj is the majority rule, 0exp and 0par were generated by means 
of an evolutionary algorithm, and 0gkl is a hand-designed rule derived from 
Gacs, Kurdyumov and Levin's work [17, 19], which is defined as following: 

majority (s*, s^.g) if s* = 
majority (s*, sj^^) ii s\ = 1 

0oxp^ uses of a block expansion [29], whereas 0par^ uses a "particle-based" 
strategy [29]. Furthermore, we considered the hand- written Das-rule [16], the 
Andre-Bennett-Koza (ABK) rule [5], which was found by means of a genetic 

2 Rule table for a r = 3 CA: 0505408305c90101200b0ef b94c7cf f 7. 
^ Rule table for a r = 3 CA: 0504058705000f 77037755837bf f b77f . 
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Table 2 



Density classification performance of CAs and RBNs for different network sizes and 
different rules. The best result in each column is shown in bold. 

algorithm, as well as two others rules generated by means of co-evolution [23]. 
Those rules will be compared to three RBN rules, namely ■d-^i'^^i'^Ti where di 
is the majority rule for K = order to compare two-states CAs with RBNs, 
the most natural way is to chose a RBN with a, K = 2r + 1 neighborhood, 
however, as we will see, even RBNs with K <2r + 1 perform better than CAs 
with a r-neighborhood. 

In [29], the performance Vn,io'^ of a network composed of nodes is defined 
as the fraction of correct classifications made over 10^ initial configurations 
that are randomly selected from an unbiased distribution, i.e., each bit in the 
initial configuration is randomly chosen. The expected percentage of Is in a 
configuration if therefore 1/2 (see Section 5). For our RBNs, VN,iQi is averaged 
over 200 randomly generated networks. 

Table 2 summarizes the performance of the density classification task for the 
different rules. One can see that the -i^-rules outperform the 0-rules since 
is as good as the best rule found by the evolutionary algorithm for A^ = 149, 
and it performs better than (/)gkl for A^ = 599 or N = 999. Furthermore, as 
mentioned in [16]: "[t]he performance of these rules decreased dramatically for 
larger A^ [•••];" which is not the case for RBNs, on the contrary, they perform 
better with increasing A^ (see Figure 4). 

In 1995, Das et al. [15] described the results obtained by four evolved 0-rules 
on the synchronization task. Those rules where used on CAs with radius r = 3 
and were evaluated by the same procedure as Mitchell et al. [29] used for the 
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Table 3 

Synchronization performance of CAs and RBNs for different network sizes and dif- 
ferent rules. 

density classification task. Table 3 summarizes the performance of 0- and 7- 
rules for different network sizes. One can see that all 7x-rules as well as the 
0sync-rule are able to perfectly solve the synchronization task. We may wonder 
why a RBN of 149 nodes does not make any error since we have seen in Figure 
10 that a network of 200 nodes did not perform perfectly? The difference lies 
in the fact that performance is essentially evaluated on configurations where 
the percentage of Is is around 1/2. 



5.4 An Entropy- Based Performance Measure 



In order to measure a rule's capability to solve the density classification or the 
synchronization task, we will introduce the following measure derived from 
Wuensche's input entropy [56]: 

N 

H{Wl I Pi) = Y.H {Wl I Pi = i/N) P (Pl = i/N) , (11) 

i=0 

where H{-) is the Shannon entropy [37]. This measure simply computes the 
conditional entropy of the number of Is in a node's input, knowing the state of 
the network at time t, i.e., the percentage of Is, Pi. Contrary to the measure 
used by Wuensche, which considers each possible node input, our measure 
only takes into account the distribution of W^, i.e., the number of Is in an 
input. This means that two node input configurations that contain the same 
number of Is are considered to be equivalent. Our measure not only allows to 
qualitatively compare the CA and RBN behavior, but it also gives indications 
on how difficult it is for a rule to solve a given task and what the average 
number of required times steps is. For example, if we obtain H{Wl \ P*) = 
after a certain amount of time, we know that the network has converged to 



22 



1.4 
1.2 
1 

0.8 
0.6 
0.4 
0.2 


2 4 6 8 10 12 14 

Fig. 12. Theoretical variation of the conditional input entropy, H{Wl \ Pf), when a 
RBN solves the density classification task. 

an all Is or all Os state. This will therefore provide us with an estimation of 
how fast a RBN may solve a task. 

We have seen in Section 3 that follows a binomial law B{K,p\), where 
p\ is implicitly defined by means of the map p\ = Qk{p^{~^) and by an initial 
distribution function for Pf. For example, Figure 12 shows the variation of the 
input entropy if we use the map defined by Equation 10 (Section 2.1) under the 
assumption that the initial configurations are randomly chosen according to a 
biased distribution, i.e., P(P{') = W(0, 1) ^ As can be seen, the input entropy 
quickly decreases towards zero. This means that — at least theoretically — it 
takes on average less than 15 time steps to settle down in a fixed point when 
K = ?) and the majority rule is used. However, as shown in Figure 5, the actual 
number of steps required is much lower (the average lies around 4.5 steps). 
Figures 13 and 14 show the variation of the conditional input entropy for the 
RBN rule "di as well as for the CA rules (^maj, 0exp, 0par, and (/>gkl {N = 149, 
initial configurations selected using a biased or unbiased distribution respec- 
tively, see Section 5.3). 

The results suggest the following comments: 

• exactly follows the curve as predicted by our theory. 

• All plots start slightly below a value of 2. This value is exactly the same 
as for a binomial law, which means that at time t = 0, no difference exists 
between RBNs and CAs. 



^ U{a, b) is the uniform distribution on the interval [a, h]. 
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Fig. 13. Variation of the conditional input entropy of different rules for the density 
classification task. Initial configurations are chosen from a biased distribution and 
= 149. (1) solid line: (j)GKL, dashed line: (/>par- 

3 I \ \ \ \ \ \ \ 1 




Fig. 14. Variation of the conditional input entropy of different rules for the density 
classification task. Initial configurations are chosen from an unbiased distribution. 
N = 149. (1) solid line: (^gkl, dashed line: (/>par- 

• When t > 0, the behaviors of the 0-rules are fairly different from the 
rules. As seen in Section 5.4, RBNs follow a simple curve that corresponds 
to a binomial law. CAs, on the other hand, need more complex behaviors 
to achieve the same result. 

• RBNs with ^^y-rule have a very simple, straightforward behavior compared 



24 



to CAs. After a very small number of iterations, the input entropy is already 
less than one bit, which means that most of the nodes basically see the same 
input. 

• The 0maj-rule is very similar to the ^^y-rule for unbiased distributions. How- 
ever, because of the CA's connection topology, frozen blocks show up and 
prevent the automata from converging to a zero input entropy. 

• The behavior of the best evolved 0par-rule and of the hand-designed (f>GKL- 
rule are very similar. Note that ^gkl has intentionally been designed to 
break the inherent symmetry of CAs and thus requires more time to con- 
verge. 

• The behavior of the 0cxp-rule is very similar to that of the ■i^y-rule. How- 
ever, Figure 14 shows that an increase in complexity is required around 
H(Wl I Pi) = 1 to overcome the emergence of frozen blocks. 0par and 0gkl, 
on the other hand, increase the complexity of their behavior from the very 
beginning in order to somehow compensate for the regularity of the inter- 
connection topology. 



Similarly, Figures 15 and 16 show the variation of the conditional input entropy 
for the synchronization task. The 77-rule as well as the three CA rules 02, 03 
and 0sync were used, = 149, and initial configurations were selected using 
a biased and unbiased distribution respectively (see Section 5.3). The results 
suggest the following comments: 



• 77 follows almost exactly the curve as predicted by our theory. A small 
difference, however, exists when the unbiased distribution is used (Figure 
16). 

• In the unbiased case, the system converges after less than three iterations, 
which is very rapid, given that the synchronization is always perfect in that 
case (see Table 3). 

• The input entropy is not fully monotonous (Figure 16). This can be ex- 
plained by looking at the evolution of the network's state during the first 
three iterations: after the first iteration, most of the nodes are in state 
since 77 produces a 1 only if the input contains Os only. The input entropy 
is thus very low, however, at the second iteration, all nodes which are con- 
nected to a 1 will be in state and all others in state 1 since they are 
connected to Os only. Thus, a single 1 in the network after the first iteration 
can generate more than one at the next step and therefore increases the 
input entropy. 

• The non-convergence of 02 and 03 after 300 iterations indicates the existence 
of frozen-blocks, which makes the problem much harder to solve for these 
rules. 
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Fig. 15. Variation of the conditional input entropy of different rules for the synchro- 
nization task. Initial configurations are chosen from a biased distribution, N = 149. 
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Fig. 16. Variation of the conditional input entropy of different rules for the syn- 
chronization task. Initial configurations are chosen from an unbiased distribution, 
N = 149. Subgraph: solid line = simulated 77, dashed line = theoretical 77. 

6 Non-Standard Architectures and Topologies 

In order to improve the computational power of CAs, various people proposed 
to introduce non-uniformity among the rules. Although this approach sacri- 
fices one of the main advantages of CAs, namely, to have simple and identical 
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processing elements, it allows to find better solutions [38,41,46]. These so- 
lutions often rely on a particular structure where the nodes' rules must be 
carefully chosen. With RBNs, however, this is not the case: as we have seen in 
Section 4, when K = b for example, a single map can represent 14, 000 rules, 
which are strictly equivalent. We can thus easily build non-uniform RBNs by 
randomly choosing one rule for each node instead of using the same rule in all 
nodes. 

CAs and RBNs can be considered as two extreme classes of networks that 
are not really observed in Nature. From the point of view of the interconnec- 
tion topology, CAs are simply too regular, whereas random networks require 
an equal probability for each connection to be established among all nodes, 
which is in most cases subjected to physical limitations since long-distance con- 
nections are more costly than short- distance connections. Many recent studies 
have instead confirmed that an important number of networks observed in Na- 
ture belong to a class called smaii- world networks (see for example [1,8,49,54]). 
These networks may be considered as an intermediate class between CAs and 
RBNs, since, if we gradually transform a CA into a RBN in a certain manner, 
we can obtain small-world networks. 

Different types of small- world networks exist [3], here we shall however only 
consider the networks of Watts and Strogatz as described in [54] . In their paper 
they state for the density classification task that "[...] a simple majority- 
rule running on a small-world graph can outperform all known human and 
genetic algorithm-generated rules running on a ring lattice." However, they 
do not give the amount of randomness a network should contain in order 
to outperform CAs on the density task. Moreover, since we also studied the 
synchronization task, it would be interesting to see if small-world networks 
using the 7-rule perform better than CAs on this task as well. We performed 
several experiments with networks built in the following way: 

(1) Start with a ring lattice, which might be considered as a special one- 
dimensional CA with a r = 1 neighborhood, but without self-connections. 

(2) Rewire the connections with probability p as described in [54]. 

Therefore, if p = 0, the initial automaton is left unchanged, if p = 1, it is trans- 
formed into an automaton with a random topology. Note that, although the 
connectivity K was 2 for each node at the beginning, each node can potentially 
have a different connectivity at the end. The average of i^T = 2 connections 
per node is, however, not affected by this random rewiring algorithm. 

Figures 17 and 18 show the percentage of success and the number of itera- 
tions required for solving the density classification and the synchronization 
task as a function of p, i.e., the amount of randomness in the network. The 
simulations were performed for a network size of TV = 149 nodes, the rules 
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Fig. 17. Percentage of success of small-world networks for the density classification 
(solid line) and the synchronization task (dashed line) as a function of the amount 
of randomness p. N = 149, majority rule (density) and 7-rule (synchronization). 



used were the majority rule for the density classification and the 7-rule for 
the synchronization task. 



The results can be summarized as following: For both tasks, the percentage 
of success reaches its maximal value around p = 0.5 (see Figure 17). Whereas 
the number of iterations for the synchronization task reaches its maximum 
around p = 0.3 (Figure 18) already, the density classification task requires 
significantly less iterations on a random topology. The reason for this can be 
found in the corresponding maps (Figures 2 and 3): the synchronization map 
does not possess stable fixed points, which makes it less sensitive to additional 
connections between the already existing clusters. 



We can conclude that a relatively small amount of randomness already greatly 
helps to improve performance for solving our two tasks. However, especially 
for the density classification task, a randomly interconnected network helps 
to significantly improve convergence. Despite the extreme simplicity of our 
random Boolean networks used, the results fit into the global picture of recent 
work on complex network synchronization using coupled oscillators (see for 
example [32,34]), although more works is certainly needed in that area. 
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Fig. 18. Number of iterations required by small-world networks for solving the den- 
sity classification (solid line) and the synchronization task (dashed line) as a func- 
tion of the amount of randomness p. N = 149, majority rule (density) and 7-rule 
(synchronization) . 

7 Conclusions 



In this paper we liave first derived a recursive equation wliicli gives the prob- 
ability of a node to be in state 1 at time t + 1 for random Boolean networks. 
Based on that equation, we were then able to analytically find the local node 
rules from the global network behavior for our two comparative tasks. Our 
main findings are the following: 

(1) No K = 2 random Boolean network can perfectly, i.e., for all initial 
configurations, solve the density classification task. The majority rule is 
the only perfect rule for K = 3 networks, whereas many rules exist for 
K > 3 networks. 

(2) The two symmetrical 7x-rules are the best rules for solving the synchro- 
nization task. 

(3) Random Boolean networks outperform cellular automata on both the 
density classification and on the synchronization task. 

(4) Random Boolean networks also outperform networks with a small-world 
topology on the density classification task, but their performance is equiv- 
alent for the synchronization task. 

(5) The rules obtained for both our tasks are highly scalable, i.e., the larger 
the network, the better they perform. This is not the case for cellular 
automata. 
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The reason why random Boolean networks perform better on both the syn- 
chronization and on the density classification tasks is the following: due to 
the random wiring, each node of a random Boolean network has — from a sta- 
tistical point of view — an unbiased view on the automata's global state since 
the distribution of Is in the node inputs of each node is the same as in the 
global network state. On the other hand, cellular automata have a biased view 
because of their local neighborhood, which makes it more difficult for them to 
solve global tasks. 

Another advantage of random Boolean networks over CAs is that they allow 
the network to be rewired at any time, without affecting the global perfor- 
mance. This property is especially interesting if interconnections are likely to 
fail. One would then only have to re-establish a new connection to a randomly 
chosen alternative node. This property could be very interesting for large scale 
distributed systems. Although not explicitly tested in this work, the existence 
of many equivalent rules for certain values of K tend to show that RBNs may 
exhibit great robustness to node and interconnection failures. 

One might of course ask whether the current approach might be applied to 
similar global tasks and whether it might be generalized to non-uniform ran- 
dom Boolean networks, alternative topologies, or asynchronously updating 
nodes. Our current work concentrates exactly on these questions. The goal 
is to extend the theory to asynchronously updating and non-uniform random 
Boolean networks and to other, more useful tasks. Asynchrony not only allows 
to remove the sole global clock signal necessary to synchronously update the 
cells, which often represents a constraint for hardware realizations, but also 
yields in further potentially interesting properties. Asynchronous automata 
attracted much interest in the past few years, although it has been proved in- 
dependently by Capcarrere [13] and Nehaniv [33] that any n-state synchronous 
automata can be emulated by a particular 3?7,^-state asynchronous automata. 
In their 1984 experimental study, Ingerson and Buvel [22] already explored the 
question of how much of the interesting behavior of cellular automata comes 
from the synchronous modeling. They concluded — based on the visual appear- 
ance of the evolution patterns over time — that the synchronous assumption 
is not essential to the study of cellular automata and that certain irrelevant 
structures may appear from the synchronous update of the cells. However, 
one of the main arguments against purely synchronous automata as a tool for 
modeling has always been the lack of biological plausibihty [48], although a 
purely asynchronous behavior is certainly not biologically plausible either. But 
there are several other issues of interest in asynchronous models: Bersini and 
Detours [9] suggested that asynchrony might induce stability, but also forces 
evolution to find more robust solutions, as Rohlfshagen and Di Paolo have 
shown for the case of rhythmic asynchronous random Boolean networks [35]. 
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A further topic of interest is also the extension of the theory to 5'-state random 
networks, which are hkely to have interesting properties and would represent a 
generalization of S'-state CAs. The formalization of asynchronously updating 
and non-uniform random Boolean or S'-state networks and the ability to derive 
their local rules for a larger class of global tasks would certainly represent a 
major step in that field of research, where the local node rules are commonly 
either hand-designed or evolved so far. 
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